#include<Algorithm/mesh_bool.h>

#include<Mesh/libcell_iterator.h>

#define quote lib_cell_quote



typedef struct My_Vertex{    
    Mesh* mesh;
    //v1 曲面mesh对应的点
    //v2 相交曲线curve对应的点
    template_v* v1,*v2;
    //store new vertex's id
    int id;
    double dis;
    void* prop;
}My_Vertex;
static inline void my_vertex_init_(My_Vertex* myv)
{
    myv->v1=NULL;myv->v2=NULL;
    myv->mesh=NULL;
    myv->id=-1;
    myv->dis=0;
    myv->prop=NULL;
}

static inline  double compute_included_angle(double *a,double*b,double*n)
{
    normalize(a,3);normalize(b,3);
    double*n1=out_product(a,b);
    double re=acos(inner_product(a,b,3));
    if(inner_product(n1,n,3)<0)
    {
        re=2*M_PI-re;
    }
    free(n1);
    return re;  
}
static inline int is_vertex_in_areas(My_Vertex**boundaries,int len,int i,My_Vertex*mv,double* normal)
{
    int j=(i-1+len)%len,k=(i+1)%len;
    if(boundaries[j]==mv||boundaries[k]==mv)
    {
        return 0;
    } 
    template_v* v=NULL;boundaries[k]->v2==NULL?(v=boundaries[k]->v1):(v=boundaries[k]->v2);
    double a[3]={v->point[0]-boundaries[i]->v2->point[0],
        v->point[1]-boundaries[i]->v2->point[1],
        v->point[2]-boundaries[i]->v2->point[2]};
    boundaries[j]->v2==NULL?(v=boundaries[j]->v1):(v=boundaries[j]->v2);
    double b[3]={v->point[0]-boundaries[i]->v2->point[0],
        v->point[1]-boundaries[i]->v2->point[1],
        v->point[2]-boundaries[i]->v2->point[2]};

    double c[3]={mv->v2->point[0]-boundaries[i]->v2->point[0],
        mv->v2->point[1]-boundaries[i]->v2->point[1],
        mv->v2->point[2]-boundaries[i]->v2->point[2]};
    double angle1=compute_included_angle(a,b,normal);
    double angle2=compute_included_angle(a,c,normal);
    if(isnan(angle1)||isnan(angle2))
    {
        return 0;
    }
    return (angle2>angle1?0:1);    
}
static inline int is_vertex_in_areas1(My_Vertex**boundaries,int len,int i,My_Vertex*mv,double* normal)
{
    printf("normal:%lf %lf %lf\n",normal[0],normal[1],normal[2]);

    int j=(i-1+len)%len,k=(i+1)%len;
    if(boundaries[j]==mv||boundaries[k]==mv)
    {
        return 0;
    } 
    template_v* v=NULL;boundaries[k]->v2==NULL?(v=boundaries[k]->v1):(v=boundaries[k]->v2);
    double a[3]={v->point[0]-boundaries[i]->v2->point[0],
        v->point[1]-boundaries[i]->v2->point[1],
        v->point[2]-boundaries[i]->v2->point[2]};
    boundaries[j]->v2==NULL?(v=boundaries[j]->v1):(v=boundaries[j]->v2);
    double b[3]={v->point[0]-boundaries[i]->v2->point[0],
        v->point[1]-boundaries[i]->v2->point[1],
        v->point[2]-boundaries[i]->v2->point[2]};

    double c[3]={mv->v2->point[0]-boundaries[i]->v2->point[0],
        mv->v2->point[1]-boundaries[i]->v2->point[1],
        mv->v2->point[2]-boundaries[i]->v2->point[2]};
    double angle1=compute_included_angle(a,b,normal);
    double angle2=compute_included_angle(a,c,normal);
    if(isnan(angle1)||isnan(angle2))
    {
        return 0;
    }
    printf("ange :%lf %lf\n",angle1,angle2);
    return (angle2>angle1?0:1);    
}


//获取cell的边界点（按顺序）
//tree1是m(c所在mesh)到新建点的映射
//tree2是相交线curve到新建点的映射

//******************
static  Node*  my_get_cell_boundary_vertices(template_c* c,Mesh2_Crossover_Point*mcp ,
    Int_RB_Tree* tree1,Int_RB_Tree*tree2)
{
    Node* re=NULL;
    // Mesh m;
    // Mesh_init(&m);
    Node* hfs=Mesh_adjust_halffaces(NULL,c);
    // Mesh_free(&m);
    for(Node* nit=hfs;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf*hf=(template_hf*)(nit->value);
        template_f* f=hf->face;
        template_v* v=hf->vertices[0];
        re=node_pushback(re,tree1->find(tree1,v->id));
        Node* node=(Node*)(mcp->f2p->find(mcp->f2p,f->id));
        for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_v*v1=(template_v*)(nit->value);
            My_Vertex* mv=(My_Vertex*)(tree2->find(tree2,v1->id));
            mv->dis=distance_of_two_points(mv->v2->point,v->point,3) ;
            for(Node* nit1=re;nit1!=NULL;nit1=(Node*)(nit1->Prev))
            {
                My_Vertex* mv1=(My_Vertex*)(nit1->value);
                if(mv1->dis<=mv->dis)
                {
                    if(re==nit1)
                    {
                        re=node_pushback(re,mv);
                    }
                    else
                    {
                        Node*nn=(Node*)malloc(sizeof(Node));
                        Node_init(nn);
                        nn->value=mv;
                        if(nit1->Next!=NULL)
                        {
                            ((Node*)(nit1->Next))->Prev=nn;
                        }
                        nn->Next=nit1->Next;
                        nit1->Next=nn;
                        nn->Prev=nit1;
                    }
                    break; 
                }
            } 
        }
    }
    free_node(hfs);
    return node_reverse(re);
}

//返回值在链表的索引
//************
static int my_get_node_value_index(Node* node,void* value)
{
    int re=-1;
    int i=0;
    for(Node*nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        if(nit->value==value)
        {
            re=i;
            break;
        }
        i++;
    }
    return re;
} 
static int array_find_value(void** a,void*v,int len)
{
    for(int i=0;i<len;i++)
    {
        if(a[i]==v)
        {
            return 1;
        }
    }

    return 0;
}
//boundaries
//normal boundaries所在平面的方向
//tree2 交线到新建点的映射(交线点到boundaries的映射)
// 给一个边界围城的区域，返回这个区域内一条切割线
static Node* my_get_split_one_split_edge(My_Vertex** boundaries,int len,double* normal,Int_RB_Tree*tree2)
{
    Node* edge=NULL;
    for(int i=0;i<len;i++)
    {
        free_node(edge);edge=NULL;
        My_Vertex* mv=boundaries[i];
        if(mv->v2==NULL)
        {
            continue;
        }
        for(Node* nit1=mv->v2->faces;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            free_node(edge);edge=NULL;
            template_f* f=(template_f*)(nit1->value);
            template_v* v=NULL;
            f->vertices[0]==mv->v2?(v=f->vertices[1]):(v=f->vertices[0]);
            My_Vertex* mv1=(My_Vertex*)(tree2->find(tree2,v->id)),*mv2=mv;
            if(mv1==NULL||is_vertex_in_areas(boundaries,len,i,mv1,normal)==0)
            {
                continue;
            }
            edge=node_pushback(edge,mv2);
            edge=node_pushback(edge,mv1);
            while(array_find_value((void**)boundaries,mv1,len)==0)
            {
                mv2=mv1;
                template_f* f1=NULL;
                for(Node*nit=mv2->v2->faces;nit!=NULL;nit=(Node*)(nit->Next))
                {
                    f1=(template_f*)(nit->value);
                    if(f1!=f)
                    {
                        f=f1;
                        break;
                    }
                } 
                if(f==NULL)
                {
                    free_node(edge);edge=NULL;
                    break;
                }
                f->vertices[0]==mv2->v2?(mv1=(My_Vertex*)(tree2->find(tree2,f->vertices[1]->id))):(mv1=(My_Vertex*)(tree2->find(tree2,f->vertices[0]->id)));
                edge=node_pushback(edge,mv1);
            }
            if(edge!=NULL)
            {
                return node_reverse(edge);
            }
        } 
    }
    return node_reverse(edge);
}

static Node* my_get_split_one_split_edge1(My_Vertex** boundaries,int len,double* normal,Int_RB_Tree*tree2)
{
    Node* edge=NULL;
    for(int i=0;i<len;i++)
    {
        free_node(edge);edge=NULL;
        My_Vertex* mv=boundaries[i];
        if(mv->v2==NULL)
        {
            continue;
        }
        //printf()
        for(Node* nit1=mv->v2->faces;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            free_node(edge);edge=NULL;
            template_f* f=(template_f*)(nit1->value);
            template_v* v=NULL;
            f->vertices[0]==mv->v2?(v=f->vertices[1]):(v=f->vertices[0]);
            My_Vertex* mv1=(My_Vertex*)(tree2->find(tree2,v->id)),*mv2=mv;
            printf("normal :%lf %lf %lf\n",normal[0],normal[1],normal[2]);
            if(mv1==NULL||is_vertex_in_areas1(boundaries,len,i,mv1,normal)==0)
            {
                continue;
            }
            printf("inere :%d %d\n",f->vertices[0]->id,f->vertices[1]->id);

            edge=node_pushback(edge,mv2);
            edge=node_pushback(edge,mv1);
            while(array_find_value((void**)boundaries,mv1,len)==0)
            {
                mv2=mv1;
                template_f* f1=NULL;
                for(Node*nit=mv2->v2->faces;nit!=NULL;nit=(Node*)(nit->Next))
                {
                    f1=(template_f*)(nit->value);
                    if(f1!=f)
                    {
                        f=f1;
                        break;
                    }
                } 
                if(f==NULL)
                {
                    free_node(edge);edge=NULL;
                    break;
                }
                f->vertices[0]==mv2->v2?(mv1=(My_Vertex*)(tree2->find(tree2,f->vertices[1]->id))):(mv1=(My_Vertex*)(tree2->find(tree2,f->vertices[0]->id)));
                edge=node_pushback(edge,mv1);
            }
            if(edge!=NULL)
            {
                return node_reverse(edge);
            }
        } 
    }
    return node_reverse(edge);
}
//给定边界围成的区域和切割线，返回切割后的两个区域
static Node** my_get_split_boundaries_area_from_boundaries_and_edge(Node* boundaries,Node* edge)
{
    if(edge==NULL)
    {return NULL;}
    Node* e_reverse=node_reverse(edge);
    int index_1=my_get_node_value_index(boundaries,edge->value);
    int index_2=my_get_node_value_index(boundaries,e_reverse->value);    
    int index1=index_1,index2=index_2;
    Node*edge1=node_copy(edge);
    if(index_1>index_2)
    {
        Node* temp_n=node_reversen(edge1);
        free_node(edge1);
        edge1=temp_n; 
        index1=index_2;
        index2=index_1;
    }
    int size=node_size(boundaries);
    Node** re=(Node**)malloc(sizeof(Node*)*2);
    memset(re,0,sizeof(Node*)*2);
    for(int i=0;i<index1;i++)
    {
        re[0]=node_pushback(re[0],node_at(boundaries,i)->value);
    }
    for(Node* nit=edge1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        re[0]=node_pushback(re[0],nit->value);
    }

    for(int i=index2+1;i<size;i++)
    {
        re[0]=node_pushback(re[0],node_at(boundaries,i)->value);
    }
    for(int i=index1+1;i<index2;i++)
    {
        re[1]=node_pushback(re[1],node_at(boundaries,i)->value);
    }
    for(Node* nit=node_reverse(edge1);nit!=NULL;nit=(Node*)(nit->Prev))
    {
        re[1]=node_pushback(re[1],nit->value);
    }
    re[0]=node_reverse(re[0]);
    re[1]=node_reverse(re[1]);
    free_node(edge1); 
    return re;
}  


//boundaries是边界
//nm是分割线
//tree2是 curve的点到boundaries点的映射
/// 给一个边界围成的区域，返回切割后的两个区域

static inline  Node** my_get_split_boundaries_areas(Node* boundaries,Mesh*curve ,Int_RB_Tree* tree2)
{
    Node** re=NULL;
    int len=node_size(boundaries);
    My_Vertex** mvs=(My_Vertex**)malloc(sizeof(My_Vertex*)*len);
    double** vv=(double**)malloc(sizeof(double*)*len);
    int i=0;
    for(Node* nit=boundaries;nit!=NULL;nit=(Node*)(nit->Next))
    {
        My_Vertex*mv=(My_Vertex*)(nit->value);
        mvs[i]=mv;
        mv->v1==NULL?(vv[i]=mv->v2->point):(vv[i]=mv->v1->point);
        i++;
    }
    double*normal=compute_plane_normal_of_jordan_curve(vv,i);
    
    Node* edge=my_get_split_one_split_edge(mvs,len,normal,tree2);
   
    re=my_get_split_boundaries_area_from_boundaries_and_edge(boundaries,edge);
    free_node(edge);
    free(mvs);free(vv);
    SAFE_FREE(normal);
    return re;
}

static inline  Node** my_get_split_boundaries_areas1(Node* boundaries,Mesh*curve ,Int_RB_Tree* tree2)
{
    Node** re=NULL;
    int len=node_size(boundaries);
    My_Vertex** mvs=(My_Vertex**)malloc(sizeof(My_Vertex*)*len);
    double** vv=(double**)malloc(sizeof(double*)*len);
    int i=0;
    for(Node* nit=boundaries;nit!=NULL;nit=(Node*)(nit->Next))
    {
        My_Vertex*mv=(My_Vertex*)(nit->value);
        mvs[i]=mv;
        mv->v1==NULL?(vv[i]=mv->v2->point):(vv[i]=mv->v1->point);
        i++;
    }
    double*normal=compute_plane_normal_of_jordan_curve(vv,i);

    printf("herer b edge\n");
    printf("normal1 :%lf %lf %lf\n",normal[0],normal[1],normal[2]);
    Node* edge=my_get_split_one_split_edge1(mvs,len,normal,tree2);
    
    printf("herer edge\n");
    re=my_get_split_boundaries_area_from_boundaries_and_edge(boundaries,edge);
    free_node(edge);
    free(mvs);free(vv);
    SAFE_FREE(normal);
    return re;
}

//curve 是相交的曲线
//mcp是m到curve的映射
//vertices储存新建的点，用来释放内存
//给一个cell,返回切割后的多边形区域
static inline int is_same_point(double *p1,double *p2)
{
    if(fabs(p1[0]-p2[0])<1e-6&&fabs(p1[1]-p2[1])<1e-6&&fabs(p1[2]-p2[2])<1e-6)
    {
        return 1;
    }
    return 0; 
}
static Node* my_get_split_areas_from_one_cell(template_c*c , Mesh2_Crossover_Point* mcp,Mesh* curve,Node** vertices)
{ 
    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    Int_RB_Tree* tree2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree2);
    for(int i=0;i<c->vertices_size;i++)
    {
        My_Vertex*mv=(My_Vertex*)malloc(sizeof(My_Vertex));
        my_vertex_init_(mv);
        *vertices=node_overlying(*vertices,mv);
        mv->v1=c->vertices[i]; 
        tree1->insert(tree1,mv->v1->id,mv);
        mv->v2=(template_v*)(mcp->v2p->find(mcp->v2p,c->vertices[i]->id));
        if(mv->v2!=NULL)
        {
            tree2->insert(tree2,mv->v2->id,mv);
            if(c->id==1132)
            {
                printf("once c v exist\n");
            }
        }   

        
    }
    for(Node* hfit=c->halffaces;hfit!=NULL;hfit=(Node*)(hfit->Next))
    {
        template_f* f=((template_hf*)(hfit->value))->face;
        Node* value_n=(Node*)(mcp->f2p->find(mcp->f2p,f->id));
        if(c->id==1132)
        {
            printf("once :%d\n",node_size(value_n));
            if(node_size(value_n)==2)
            {
                if(is_same_point( ((template_v*)(value_n->value) )->point,
                    ((template_v*)( ((Node*)(value_n->Next))->value ) )->point  ) )
                {
                    printf("is_same_point\n");
                }
                else
                {
                    template_v* v1=((template_v*)(value_n->value) ),
                        *v2=((template_v*)( ((Node*)(value_n->Next))->value ) );
                    printf("not_same_point %d %d\n",v1->id,
                        v2->id);
                    printf("%d %d\n",node_size(v1->faces),node_size(v2->faces) );
                    for(Node* nit=v1->faces;nit!=NULL;nit=(Node*)(nit->Next))
                    {
                        f=(template_f*)(nit->value);
                        printf("%d %d\n",f->vertices[0]->id,f->vertices[1]->id );

                    }
                    for(Node* nit=v2->faces;nit!=NULL;nit=(Node*)(nit->Next))
                    {
                        f=(template_f*)(nit->value);
                        printf("%d %d\n",f->vertices[0]->id,f->vertices[1]->id );

                    }
                }
            }

        }
        for(Node* nit=value_n;nit!=NULL;nit=(Node*)(nit->Next))
        {
            My_Vertex* mv=(My_Vertex*)malloc(sizeof(My_Vertex));
            my_vertex_init_(mv);
            *vertices=node_overlying(*vertices,mv);
            mv->v2=(template_v*)(nit->value);

            if(c->id==1132)
            {
                printf("id:%d %lf %lf %lf\n",mv->v2->id,mv->v2->point[0],mv->v2->point[1],mv->v2->point[2]  );
            }

            tree2->insert(tree2,mv->v2->id,mv);
        }
    }

    Node* node=(Node*)(mcp->c2p->find(mcp->c2p,c->id));
    if(c->id==1132)
    {
        printf("once c  exist %d\n",node_size(node));
    }
    for(Node* nit=node; nit!=NULL;nit=(Node*)(nit->Next))
    {

        My_Vertex* mv=(My_Vertex*)malloc(sizeof(My_Vertex));
        my_vertex_init_(mv);
        *vertices=node_overlying(*vertices,mv);
        mv->v2=(template_v*)(nit->value);
        if(c->id==1132)
        {
            printf("id:%d %lf %lf %lf\n",mv->v2->id,mv->v2->point[0],mv->v2->point[1],mv->v2->point[2]  );
        }
        tree2->insert(tree2,mv->v2->id,mv); 
    } 
    if(c->id==1132)
    {
        printf("split hrere: %d %d %d\n",c->vertices[0]->id,c->vertices[1]->id,c->vertices[2]->id );
        printf("size:%d\n",node_size(*vertices));
    }
    Node* re=NULL,*re1=NULL;

    Node* boundaries=my_get_cell_boundary_vertices(c,mcp,tree1,tree2);
    if(c->id==1132)
    {
        printf("split hrere1 %d \n",node_size(boundaries));

        for(Node* nit=boundaries;nit!=NULL;nit=(Node*)(nit->Next))
        {
            My_Vertex* mv=(My_Vertex*)(nit->value);
            if(mv->v1!=NULL)
            {
                printf("v1;%d\n",mv->v1->id);
            }
            else
            {
                printf("v2;%d\n",mv->v2->id);
            }

        }
    }
    re1=node_overlying(re1,boundaries);
    while(re1!=NULL)
    { 
        Node* temp_n=NULL;
        for(Node* nit=re1;nit!=NULL;nit=(Node*)(nit->Next))
        { 
            Node**two_n=NULL;
            if(c->id==1132)
            {
                printf("begin two\n");


                for(Node* nit1=(Node*)(nit->value);nit1!=NULL;nit1=(Node*)(nit1->Next))
                {
                    My_Vertex* mv=(My_Vertex*)(nit1->value);
                    if(mv->v1!=NULL)
                    {
                        printf("v1;%d\n",mv->v1->id);
                    }
                    else
                    {
                        printf("v2;%d\n",mv->v2->id);
                    }

                }
                two_n= my_get_split_boundaries_areas1((Node*)(nit->value),curve,tree2);

            }
            else
            {
                two_n= my_get_split_boundaries_areas((Node*)(nit->value),curve,tree2);
            }

            if(c->id==1132)
            {
                printf("end wto\n");

            }

            if(two_n==NULL)
            {
                re=node_overlying(re,nit->value);
            }
            else
            {   
                temp_n=node_overlying(temp_n,two_n[0]);
                temp_n=node_overlying(temp_n,two_n[1]);  
                free(two_n);
                free_node((Node*)(nit->value));
            }
        }
        free_node(re1);
        re1=temp_n;
    } 
    if(c->id==1132)
    {
        printf("split hrere2\n");
    }
    int_rb_tree_free(tree1);int_rb_tree_free(tree2);
    return re;
} 


//如果稳定就修改,主要对那些没有交线的cell直接返回。

Mesh* remesh_along_intersecting_curve(Node* m1,Mesh2_Crossover_Point*mcp,Mesh*curve,Int_RB_Tree* tree2,Int_RB_Tree* tree3)
{
    if(m1==NULL)
    {
        return NULL;
    } 
    Mesh* re=(Mesh*)malloc(sizeof(Mesh));
    Mesh_init(re);
    re->dimension=2;re->simplex=1;re->manifold_require=1;
    //tree1是m1的点到re的点映射

    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    for(Node* cit=m1;cit!=NULL;cit=(Node*)(cit->Next))
    { 
        Node* vertices=NULL;
        //printf("begin spllit area \n");
        printf("remsh c id:%d\n",((template_c*)(cit->value))->id);
        
        Node* bns=my_get_split_areas_from_one_cell((template_c*)(cit->value),mcp,curve,&vertices); 
        

        if( ((template_c*)(cit->value))->id==1132)
        {
            printf("split areas :%d\n",node_size(bns));
        } 
        for(Node* nit=bns;nit!=NULL;nit=(Node*)(nit->Next))
        { 
            Node* node=(Node*)(nit->value);
            int len=node_size(node); 
            if( ((template_c*)(cit->value))->id==1132)
            {
                printf("sares size:%d\n",len);
                //printf("split areas :%d\n",node_size(bns));
            } 
            double**vv=(double**)malloc(sizeof(double*)*len);
            Node* nit1=node;
            template_v* v=NULL; 
            for(int i=0;i<len;i++)
            {
                vv[i]=(double*)malloc(sizeof(double)*3);
                My_Vertex* mv=(My_Vertex*)(nit1->value);
                v=(mv->v2==NULL?mv->v1:mv->v2); 

                if( ((template_c*)(cit->value))->id==1132)
                {
                    if(mv->v1!=NULL)
                    {
                        printf("v1 %d %.9lf %.9lf %.9lf ",mv->v1->id,v->point[0],v->point[1],v->point[2]);
                    }
                    else
                    {
                        printf("v2 %d %.9lf %.9lf %.9lf ",mv->v2->id,v->point[0],v->point[1],v->point[2]);
                    }
                //printf("split areas :%d\n",node_size(bns));
                } 
                vv[i][0]=v->point[0];vv[i][1]=v->point[1];vv[i][2]=v->point[2]; 
                nit1=(Node*)(nit1->Next);
            }
            int** s=( int**)malloc(sizeof(int*)*(len-2));
            for(int i=0;i<len-2;i++)
            {
                s[i]=(int*)malloc(sizeof(int)*3);
                memset(s[i],0,sizeof(int)*3);
            }
            optimization_subdivision_of_polygon(vv,len,s);
            if( ((template_c*)(cit->value))->id==1132)
            {
                printf("begin optimizat\n");
                //printf("split areas :%d\n",node_size(bns));
            } 
            for(int i=0;i<len-2;i++)
            { 
                template_v* vs[3]={NULL,NULL,NULL};
                for(int j=0;j<3;j++)    
                {
                    My_Vertex* mv=(My_Vertex*)(node_at(node,s[i][j])->value);
                    vs[j]=(mv->v2==NULL?((template_v*)(tree1->find(tree1,mv->v1->id))):((template_v*)(tree2->find(tree2,mv->v2->id))));
                    if(vs[j]==NULL)
                    { 
                        vs[j]=re->create_vertexv(re,(mv->v2==NULL?mv->v1->point:mv->v2->point),3);
                        mv->v2==NULL?tree1->insert(tree1,mv->v1->id,vs[j]):tree2->insert(tree2,mv->v2->id,vs[j]); 
                    }
                    if( ((template_c*)(cit->value))->id==1132)
                    {
                        if(mv->v1!=NULL)
                        {
                            printf("v1 %d %d",mv->v1->id,s[i][j]);
                        }
                        else
                        {
                            printf("v2 %d %d",mv->v2->id,s[i][j]);
                        }
                    } 
                   
                }
                if( ((template_c*)(cit->value))->id==1132)
                {
                    printf("\n");
                    if(vs[0]==NULL||vs[1]==NULL||vs[2]==NULL)
                    {
                        printf("cuowu\n");
                    }
                //printf("split areas :%d\n",node_size(bns));
                } 

                template_c*cc=re->create_cellv(re,vs,3); 
                if( ((template_c*)(cit->value))->id==1132)
                {
                    printf("end create cell\n");
                //printf("split areas :%d\n",node_size(bns));
                } 

                if(cc!=NULL)
                {
                    tree3->insert(tree3,cc->id,cit->value);
                }
                else
                {
                    // template_v*tempvvv=vs[0];
                    // vs[0]=vs[1];vs[1]=tempvvv;
                    // cc=re->create_cellv(re,vs,3); 
                    // if(cc!=NULL)
                    // {
                    //     tree3->insert(tree3,cc->id,cit->value);

                    // }

                    // if(cc==NULL)
                    // {
                        printf("zheli c id:%d 错 误  return NULL %d \n",((template_c*)(cit->value))->id,len);
                        //Mesh_free(re);free(re);

                    //printf("zheli c修改\n");
                    //     return NULL;
                    // }
                   
                }
                if( ((template_c*)(cit->value))->id==1132)
                {
                    printf("end hrere\n");

                }

            } 
            if( ((template_c*)(cit->value))->id==1132)
            {
                printf("end hrere1\n");

            }

            free_node(node);
            for(int i=0;i<len;i++)
            {
                free(vv[i]);
            }
            for(int i=0;i<len-2;i++)
            {
                free(s[i]);
            }
            free(vv);free(s);
            if( ((template_c*)(cit->value))->id==1132)
            {
                printf("end hrere2\n");
            }

        } 
        if( ((template_c*)(cit->value))->id==1132)
        {
            printf("hree4\n");
        } 
        free_node_value(vertices);
        free_node(vertices);
        free_node(bns);
        if( ((template_c*)(cit->value))->id==1132)
        {
            printf("hree5\n");
        } 
    }
    int_rb_tree_free(tree1); 
    return re;
}

//curve 是相交的曲线
//mcp是m到curve的映射
//vertices储存新建的点，用来释放内存
//给一个cell,返回切割后的多边形区域


static Node* my_get_sub_info_from_one_cell(template_c* c, Mesh2_Crossover_Point* mcp, Mesh* curve, 
    Node** vertices, Node** inner_edges)
{
    int id = 0;
    Int_RB_Tree* tree1 = (Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    Int_RB_Tree* tree2 = (Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree2);
    for (int i = 0; i < c->vertices_size; i++)
    {
        My_Vertex* mv = (My_Vertex*)malloc(sizeof(My_Vertex));
        my_vertex_init_(mv); mv->id = id; id++;
        *vertices = node_overlying(*vertices, mv);
        mv->v1 = c->vertices[i];
        tree1->insert(tree1, mv->v1->id, mv);
        mv->v2 = (template_v*)(mcp->v2p->find(mcp->v2p, c->vertices[i]->id));
        if (mv->v2 != NULL)
        {
            tree2->insert(tree2, mv->v2->id, mv);
        }
    }
    Int_RB_Tree* tree3 = (Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree3);
    Int_RB_Tree* tree4 = (Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree4);
    //printf("begin  sub info\n");
    for (Node* hfit = c->halffaces; hfit != NULL; hfit = (Node*)(hfit->Next))
    {
        template_f* f = ((template_hf*)(hfit->value))->face;
        Node* value_n = (Node*)(mcp->f2p->find(mcp->f2p, f->id));
        //printf("f id :%d\n",f->id);
        for (Node* nit = value_n; nit != NULL; nit = (Node*)(nit->Next))
        {
            My_Vertex* mv = (My_Vertex*)malloc(sizeof(My_Vertex));
            my_vertex_init_(mv); mv->id = id; id++;
            *vertices = node_overlying(*vertices, mv);

            mv->v2 = (template_v*)(nit->value);
            tree4->insert(tree4,mv->v2->id,mv->v2);
            //printf("insert tree4 id :%d\n",mv->v2->id);
            for(Node* nit1=mv->v2->faces;nit1!=NULL;nit1=(Node*)(nit1->Next))
            {
                //*********
                template_f* f1=(template_f*)(nit1->value);
                //printf("f1 %d %d\n",f1->vertices[0]->id,f1->vertices[1]->id);
                template_v* v=NULL;
                f1->vertices[0]==mv->v2?(v=f1->vertices[1]):(v=f1->vertices[0]);
                Crossover_Point* cp=(Crossover_Point*)(v->prop);
                if(cp->scp[0].f!=f&&cp->scp[1].f!=f&&tree4->find(tree4,v->id)!=NULL)
                {
                    //printf("insert\n");
                    tree3->insert(tree3,f1->id,f1);
                }
            }
            tree2->insert(tree2, mv->v2->id, mv);
        }
    }

    Node* node = (Node*)(mcp->c2p->find(mcp->c2p, c->id));
   
    for (Node* nit = node; nit != NULL; nit = (Node*)(nit->Next))
    {
        My_Vertex* mv = (My_Vertex*)malloc(sizeof(My_Vertex));
        my_vertex_init_(mv); mv->id = id; id++;
        *vertices = node_overlying(*vertices, mv);
        mv->v2 = (template_v*)(nit->value);
      
        tree2->insert(tree2, mv->v2->id, mv);
    }

    

    for (Node* nit = node; nit != NULL; nit = (Node*)(nit->Next))
    {
        for (Node* nit1 = ((template_v*)(nit->value))->faces; nit1 != NULL; nit1 = (Node*)(nit1->Next))
        {
            template_f* f = (template_f*)(nit1->value);
            tree3->insert(tree3, f->id, f);
            
        }
    }
    for (auto it=tree3->begin(tree3);it.it!=NULL;it++)
    {
        template_f* f = (template_f*)(it.second);
        //tree3->insert(tree3, f->id, f);
        My_Vertex** temp = (My_Vertex**)malloc(sizeof(My_Vertex*) * 2);
        memset(temp, 0, sizeof(My_Vertex*) * 2);
        for (int i = 0; i < 2; i++)
        {
            temp[i] = (My_Vertex*)(tree2->find(tree2, f->vertices[i]->id));
            if (temp[i] == NULL)
            {
                temp[i] = (My_Vertex*)(tree1->find(tree1, f->vertices[i]->id));
            }
            if(temp[i]==NULL)
            {
                SAFE_FREE(temp); break;
            }
        }
        if(temp!=NULL)
        {
            *inner_edges = node_overlying(*inner_edges, temp);
        }
    }

    Node* boundaries = my_get_cell_boundary_vertices(c, mcp, tree1, tree2);
     
    int_rb_tree_free(tree1); int_rb_tree_free(tree2);
    int_rb_tree_free(tree3);    int_rb_tree_free(tree4);

    return boundaries;
}


//m1是要remesh的网格区域
//curve 
//mcp是m1到curve的映射
//tree2 是curve的点到返回mesh的点映射
//tree3 是返回mesh的cell到 m1的cell映射。 

Mesh* remesh_along_intersecting_curven(Node* m1, Mesh2_Crossover_Point* mcp, Mesh* curve, Int_RB_Tree* tree2, Int_RB_Tree* tree3)
{
    if (m1 == NULL)
    {
        return NULL;
    }
    
    Mesh* re = (Mesh*)malloc(sizeof(Mesh));
    Mesh_init(re);
    re->dimension = 2; re->simplex = 1; re->manifold_require = 1;
    //tree1是m1的点到re的点映射

    Int_RB_Tree* tree1 = (Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    //printf("begin remesh 1\n");
    for (Node* cit = m1; cit != NULL; cit = (Node*)(cit->Next))
    {
        Node* vertices = NULL, *innder_edges = NULL;
        //printf("begin spllit area \n");
        //printf("remsh c id:%d\n", ((template_c*)(cit->value))->id);
        // if (((template_c*)(cit->value))->id != 7378)
        // {
        //     continue;
        // }
        if(((template_c*)(cit->value))->id==7378)
        {
            printf("remsh c id:%d\n", ((template_c*)(cit->value))->id);
        }
        //printf("remsh c id:%d\n", ((template_c*)(cit->value))->id);
        Node* boundaries = my_get_sub_info_from_one_cell((template_c*)(cit->value), mcp, curve, &vertices,&innder_edges);
        //printf("once end boundaries\n");
        int vs_len = node_size(vertices);
        int len_bounds = node_size(boundaries);
        int inner_len = node_size(innder_edges);
        double** vs = (double**)malloc(sizeof(double*)*vs_len);
        My_Vertex** vs1 = (My_Vertex**)malloc(sizeof(My_Vertex*)*vs_len);
        for (Node* nit=vertices;nit!=NULL;nit=(Node*)(nit->Next))
        {
            My_Vertex* mv = (My_Vertex*)(nit->value);
            vs1[mv->id] = mv;
            mv->v2 == NULL ? (vs[mv->id] = mv->v1->point) : (vs[mv->id] = mv->v2->point);

            if (((template_c*)(cit->value))->id == 7378)
            {
                if (mv->v2!=NULL)
                {
                    printf("v2 :%d %d %lf %lf %lf\n", mv->v2->id,mv->id,
                        mv->v2->point[0],mv->v2->point[1],mv->v2->point[2]);
                }
                else
                {
                    printf("v1 :%d %d %lf %lf %lf\n", mv->v1->id,mv->id,
                        mv->v1->point[0],mv->v1->point[1],mv->v1->point[2]);
                }
            }
        }
        int* bounds = NULL;
        int** inner_edges1 = NULL;

        if (len_bounds > 0)
        {
            if (((template_c*)(cit->value))->id == 7378)
            {
                printf("len boudns :%d\n",len_bounds);
            }
            int i = 0;
            bounds = (int*)malloc(sizeof(int)*len_bounds);
            for (Node* nit= boundaries;nit!=NULL;nit=(Node*)(nit->Next))
            {
                My_Vertex* mv = (My_Vertex*)(nit->value);
                if (((template_c*)(cit->value))->id == 7378)
                {
                    if (mv->v2!=NULL)
                    {
                        printf("v2 :%d %d %lf %lf %lf\n", mv->v2->id,mv->id,
                        mv->v2->point[0],mv->v2->point[1],mv->v2->point[2]);
                    }
                    else
                    {
                        printf("v1 :%d %d %lf %lf %lf\n", mv->v1->id,mv->id,
                        mv->v1->point[0],mv->v1->point[1],mv->v1->point[2]);
                    }
                }
                bounds[i] = mv->id;
                i++;
            }
        }
        if (inner_len > 0)
        {
            if (((template_c*)(cit->value))->id == 7378)
            {
                printf("inner len :%d\n",inner_len);
            }
            int i = 0;
            inner_edges1 = (int**)malloc(sizeof(int*) * inner_len);
            for (Node* nit = innder_edges; nit != NULL; nit = (Node*)(nit->Next))
            {
                My_Vertex** mv = (My_Vertex**)(nit->value);
                inner_edges1[i] = (int*)malloc(sizeof(int) * 2);
                inner_edges1[i][0] = mv[0]->id;inner_edges1[i][1] = mv[1]->id;
                if (((template_c*)(cit->value))->id == 7378)
                {
                    printf("edge %d %d\n",mv[0]->id,mv[1]->id);
                }
                i++;
            }
        }
        // if (((template_c*)(cit->value))->id == 5351)
        // {
        //     printf("inner_len: %d\n",inner_len);
        // }
        //printf("sub d begin\n");
        Node*tris= subdivision_of_polygon_with_edges1(vs, vs_len,  bounds, len_bounds,
            inner_edges1, inner_len);
        //printf("sub d end\n" );
        //int flag=0;
        for (Node* nit=tris;nit!=NULL;nit=(Node*)(nit->Next))
        {
            int* tri = (int*)(nit->value);
            if (((template_c*)(cit->value))->id == 7378)
            {
                printf("tri:%d %d %d\n",tri[0],tri[1],tri[2] );
            }
            template_v* vvs[3] = {NULL,NULL,NULL};
            for (int j=0;j<3;j++)
            {
                My_Vertex* mv = (My_Vertex*)(vs1[tri[j]]);
                vvs[j] = (mv->v2 == NULL ? ((template_v*)(tree1->find(tree1, mv->v1->id))) : ((template_v*)(tree2->find(tree2, mv->v2->id))));
                
                if (vvs[j] == NULL)
                {
                    vvs[j] = re->create_vertexv(re, (mv->v2 == NULL ? mv->v1->point : mv->v2->point), 3);
                    mv->v2 == NULL ? tree1->insert(tree1, mv->v1->id, vvs[j]) : tree2->insert(tree2, mv->v2->id, vvs[j]);
                    
                    // if (mv->v2 != NULL&&(mv->v2->id==47||mv->v2->id==43) )
                    // {
                    //     //flag=1;
                    //     printf("inset v2:%d %d\n ", mv->v2->id,vvs[j]->id);
                    //     printf("remsh c id:%d\n", ((template_c*)(cit->value))->id);
 
                    // }
                    // else
                    // {
                    //     printf("inset v1:%d\n ", mv->v1->id);
                    // }
                }
                // if (mv->v2 != NULL&&(mv->v2->id==47||mv->v2->id==43) )
                // {
                //         //flag=1;
                //     printf("inset v2:%d %d\n ", mv->v2->id,vvs[j]->id);
                //     printf("remsh c id:%d\n", ((template_c*)(cit->value))->id);
                // }
            }

            template_c* cc = re->create_cellv(re, vvs, 3);
            // if(flag)
            // {
            //     printf("%d %d %d c :%d\n",vvs[0]->id,vvs[1]->id,vvs[2]->id,((template_c*)(cit->value))->id);
            // }
            if (cc != NULL)
            {
                //printf("inber\n");
                tree3->insert(tree3,cc->id,cit->value);
            }
            else
            {
                printf("cuowu ********************* %d\n", ((template_c*)(cit->value))->id);
                return re;
            }

        }

        free(vs1); free(vs);
        free_node_value(tris);
        free_node(tris);
        SAFE_FREE(bounds);
        for (int i=0;i<inner_len;i++)
        {
            free(inner_edges1[i]);
        }
        SAFE_FREE(inner_edges1);
        free_node_value(vertices);
        free_node(vertices);
        free_node(boundaries);
        free_node_value(innder_edges);
        free_node(innder_edges);

    }
    int_rb_tree_free(tree1);
    return re;
}

static Node* get_union_faces_boundaries(Node* node1)
{
    Node* re=NULL;

    while(node1!=NULL)
    {
        Node* node=NULL;
        template_f* f=(template_f*)(node1->value);
        node1=(Node*)(node1->Next);
        node=node_overlying(node,f->vertices[0]);
        node=node_overlying(node,f->vertices[1]);
        int flag=1;
        while(flag==1)
        {
            flag=0;
            template_v* v1=(template_v*)(node->value);
            template_v* v2=(template_v*)(node_reverse(node)->value);
            for(Node* nit=node1;nit!=NULL;nit=(Node*)(nit->Next))
            {
                template_f* f1=(template_f*)(nit->value);
                if(f1->vertices[0]==v1)
                {
                    node=node_overlying(node,f1->vertices[1]);
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
                else if(f1->vertices[1]==v1)
                {
                    node=node_overlying(node,f1->vertices[0]);
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
                else if(f1->vertices[0]==v2)
                {
                    node_pushback(node_reverse(node),f1->vertices[1]);
                    
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
                else if(f1->vertices[1]==v2)
                {
                    node_pushback(node_reverse(node),f1->vertices[0]);
                    
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
            }
        }
        re=node_overlying(re,node);
    }


    return re;
}

void cut_along_curve_mesh(Mesh* m,Mesh* curve,Int_RB_Tree* tree,Int_RB_Tree* tree0,Int_RB_Tree* tree4,  Int_RB_Tree *tree3)
{
    if(m==NULL)
    {
        return;
    }

    Node* node=NULL,*node1=NULL;

    for(auto fit=curve->f_begin(curve);fit!=curve->f_end(curve);fit++)
    {
        template_v* vs[2]={NULL};
        vs[0]=(template_v*)(tree->find(tree,quote(fit)->vertices[0]->id));
        vs[1]=(template_v*)(tree->find(tree,quote(fit)->vertices[1]->id));
        if(vs[0]==NULL||vs[1]==NULL)
        {
            node1=node_overlying(node1,quote(fit));
            printf("zheli f 修改\n");
            if(vs[0]==NULL)
            {
                printf("v id:%d\n",quote(fit)->vertices[0]->id);
            } 
            if(vs[1]==NULL)
            {
                printf("v id:%d\n",quote(fit)->vertices[1]->id);
            } 
            continue;
        }
        template_f* f=m->get_facev(m,vs,2);
        if(f==NULL)
        {

            printf("cuowu vs:%d %d from :%d %d\n",vs[0]->id,vs[1]->id,
                quote(fit)->vertices[0]->id,quote(fit)->vertices[1]->id);

            Node*fvs= get_shortest_path(m, vs[0],vs[1]);
            for(Node* nit=fvs;nit!=NULL;nit=(Node*)(nit->Next))
            {
                if(nit->Next!=NULL)
                {
                    vs[0]=(template_v*)(nit->value);
                    vs[1]=(template_v*)((Node*)(( (Node*)(nit->Next) )->value));
                    node=node_overlying(node,m->get_facev(m,vs,2) );
                }

                //printf("%d\n",((template_v*)(nit->value))->id );
            }
            free_node(fvs);

        }
        else
        {
            node=node_overlying(node,f);
        }

    }  
    Node* unios_bounds=get_union_faces_boundaries(node1);
    for(Node* nit=unios_bounds;nit!=NULL;nit=(Node*)(nit->Next))
    {   
        Node* node2=(Node*)(nit->value);
        template_v* vs[2]={0};
        vs[0]= (template_v*)(tree->find(tree,((template_v*)(node2->value))->id)) ;
        vs[1]=(template_v*)(tree->find(tree,((template_v*)(node_reverse(node2)->value))->id));
        //printf("unios :%d :%d\n",vs[0]->id,vs[1]->id);
        Node*fvs= get_shortest_path(m, vs[0],vs[1]);
        for(Node* nit1=fvs;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            if(nit1->Next!=NULL)
            {
                vs[0]=(template_v*)(nit1->value);
                vs[1]=(template_v*)((Node*)(( (Node*)(nit1->Next) )->value));
                node=node_overlying(node,m->get_facev(m,vs,2) );
            }

           // printf("%d\n",((template_v*)(nit->value))->id );
        }
        free_node(fvs);

        free_node(node2);
    }
    free_node(unios_bounds);
    //printf("hereend\n");
    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    Int_RB_Tree* tree2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree2);
    Mesh_cut_along_the_curvef(m,node,tree1,tree2);

    if(tree4!=NULL)
    {
        for(auto it=tree2->begin(tree2);it.it!=NULL;it++)
        {
            template_v* v=(template_v*)(it.second);
            template_v* vv=(template_v*)(tree0->find(tree0,v->id ));
            if(vv==NULL)
            {
                printf("zheli vv 修改\n");
            }
            else
            {
                tree4->insert(tree4,it.first,vv);
            }
        }
    }
    
    //printf("end cut cuvef\n");
    if(tree3!=NULL)
    {
        for(auto it=tree1->begin(tree1);it.it!=NULL;it++)
        {
            int id=*((int*)(it.second));

            void* value=tree3->find(tree3,id);
            if(value!=NULL)
            {
                tree3->erase(tree3,id);
                tree3->insert(tree3,it.first,value); 
            }
        }
    }
    int_rb_tree_free_value(tree1);
    int_rb_tree_free(tree1); 
    int_rb_tree_free(tree2);
    //int_rb_tree_free(tree0);
    free_node(node);
}


